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We present the first results from our new general relativistic, Lagrangian hydrodynamics code, 
which treats gravity in the conformally flat (CF) limit. The evolution of fluid configurations is 
described using smoothed particle hydrodynamics (SPH), and the elliptic field equations of the 
CF formalism are solved using spectral methods in spherical coordinates. The code was tested on 
models for which the CF limit is exact, finding good agreement with the classical Oppenheimer- 
Volkov solution for a relativistic static spherical star as well as the exact semi-analytic solution 
for a collapsing spherical dust cloud. By computing the evolution of quasi-equilibrium neutron 
star binary configurations in the absence of gravitational radiation backreaction, we have confirmed 
that these configurations can remain dynamically stable all the way to the development of a cusp. 
With an approximate treatment of radiation reaction, we have calculated the complete merger of 
an irrotational binary configuration from the innermost point on an equilibrium sequence through 
merger and remnant formation and ringdown, finding good agreement with previous relativistic 
calculations. In particular, we find that mass loss is highly suppressed by relativistic effects, but 
that, for a reasonably stiff neutron star equation of state, the remnant is initially stable against 
gravitational collapse because of its strong differential rotation. The gravity wave signal derived from 
our numerical calculation has an energy spectrum which matches extremely well with estimates based 
solely on quasi-equilibrium results, deviating from the Newtonian power-law form at frequencies 
below 1 kHz, i.e., within the reach of advanced interferometric detectors. 



I. INTRODUCTION 



Gravitational wave (GW) astronomy stands at a crucial moment in its history, with the LIGO (Laser Interferometer 
Gravitational Wave Observatory) Scientific Collaboration reporting results from their first scientific runs Q, 0, El , 
GEO600 having completed two scientific runs 0, Q , TAMA taking data 0, , and VIRGO reporting its first lock 
acquisition 0, U ITol ] . As such, it is now more important than ever to have accurate theoretical predictions of the 
main candidate signals, both to aid in their detection and to facilitate the interpretation of any future detections. 

It has long been recognized that coalescing relativistic binary systems containing compact objects, either neutron 
stars (NS) or black holes (BH), are likely to be important sources of detectable GWs. Recent population synthesis 
calculations indicate that an Advanced LIGO detector should be able to see at least tens of coalescences per year of 
NS-NS, NS-BH, and BH-BH binaries 11]. Empirical rate estimates based on three of the four observed binary pulsar 
systems expected to merge within a Hubble time, PSR B1913+16, PSR B1534+12, and PSR J0737-3039 pjj, are in 
general agreement, with the latter making the dominant contribution to the probabilistic rate because of its short 
coalescence time of 85 Myr (see 0, 0] and references therein; the fourth system is located in a globular cluster, 
rather than the galactic plane). NS-NS binaries are the only known systems with coalescence times shorter than a 
Hubble time to have been conclusively observed, and it is the orbital decay of PSR1913+16 that currently provides 
our best indirect evidence for the existence of GWs [TEl [Til Il7j . 

At large separations, the dynamics of compact object binaries can be well approximated by high-order post- 
Newtonian (PN) (see ^^] and references therein) and other formalisms |2(j which treat the compact objects 
as point-masses, including the effects of spin-orbit and spin-spin angular momentum couplings for systems containing 
a BH [2ll|23,|23. In general, these methods are more appropriate for describing a BH than a NS, since the finite-size 
effects ignored by point-mass formulae are of greater magnitude in systems containing NS. 

For smaller binary separations, r < IORns, where Rns is the radius of the NS, these finite-size corrections play an 
increasingly significant role in the evolution of NS-NS binaries. As long as the coalescence timescale r = r/f remains 
longer than the dynamical timescale, the evolution is quasi-adiabatic, and the binary will sweep through a sequence 
of configurations representing energy minima for given binary separations. The infall rate will be given by 
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where (dE/dt)aw is the energy loss rate to gravitational radiation, and (dE/dr) eq is the slope of the equilibrium 
energy curve. Equilibrium energy sequences were first constructed in Newtonian and then PN gravity (see |2 1| and 
references therein). More recently, general relativistic (GR) sequences have been calculated for binary NS systems 
in quasi-circular orbits l2(jL WA Esl I29I I30I l3ll l32] | . An important result from these studies is that the slope of 
the equilibrium energy curve is flatter when relativistic effects are included, compared to the Newtonian value. This 
becomes especially pronounced at small separations, near the point where equilibrium sequences encounter an energy 
minimum (often associated with the "innermost stable circular orbit" , or ISCO) or terminate when a cusp develops on 
the inner edge of the NS. Based on this observation, it was noted that a flatter slope in the equilibrium energy curve 
results in not only a faster infall rate but also a decrease in the energy spectrum dE/df ~ (dE / dr) / (df / dr) . In fact, 
based on the relativistic equilibrium sequences of Taniguchi and Gourgoulhon |27l hereafter TG], Faber et al. |33l 
hereafter FGRT] concluded that the dependence of finite-size corrections on the NS compactness (i.e. M^g/R^s) 
would leave an imprint in the GW energy spectrum at frequencies / « 500 — 1000 Hz, within the band accessible to 
Advanced LIGO. Hughes [34[ proposed a method using these results which would allow for a determination of the NS 
compactness to within a few percent based on 10 — 50 observations with Advanced LIGO if it employs a narrow-band 
detector in addition to the current broadband setup, with the required number of observations dependent upon both 
the true compactness and the particular setup of the detector. 

Shortly before the ISCO or the end of the equilibrium sequence (when it terminates at the formation of a cusp) the 
binary will begin a transition toward a rapid plunge inward, eventually leading to merger. Once the binary passes 
this point, the dynamical evolution becomes too complicated to describe using semi-analytical methods, requiring 
instead a 3-d hydrodynamic treatment. Such calculations were performed first for Newtonian gravitation, using both 
Eulerian, grid-based codes and Lagrangian smoothed particle hydrodynamics (SPH; for a review, see, e.g., |35j | and 
references therein). It was always recognized that any results would be at best qualitatively accurate, since the extreme 
compactness of the NS would induce a host of GR effects. No ting this, some attempts were made to calculate the 
evolution of binary NS systems in lowest-order PN gravity [3^, l36t l37l l38l l39| , using a formalism developed by 0] • 
Unfortunately, using realistic NS equations of state (EOS) violated the basic assumption of the PN approximation 
that the magnitude of the 1PN terms are small relative to Newtonian-order effects. As a result, all PN calculations 
were forced to make unphysical approximations, either by evolving NS with a fraction of their proper physical mass 
[33,133, or by reducing the magnitude of all 1PN terms (3^, [S^, |33> hereafter denoted FR1-3, or collectively FR]. 

While the ultimate goal of studying binary NS coalescences should be a full GR treatment, only one group has 
been able to calculate the full evolution of a binary system from an equilibrium binary configuration through merger 
and the formation of a remnant [4ll I42I Elij . While these calculations represent a breakthrough in our understanding 
of the hydrodynamics of coalescing NS binaries, they still leave a great deal of room for further research into a 
variety of questions. Currently, full GR calculations are extremely difficult, and are vulnerable to several numerical 
instabilities. In order to guarantee the stability of calculation through coalescence and the formation of either a 
stable merger remnant or a BH, binaries were started from the termination points of equilibrium sequences, from 
quasi-circular configurations with zero infall velocity. The calculations of Shibata and Uryu [ill |42| used initial 
conditions generated by Uryu, Eriguchi, and collaborators |3lll32|. and the more recent work by Shibata, Taniguchi, 
and Uryu |43j used initial conditions generated by TG. This restricts studying the behavior of the infall velocity and 
GW signal when the stars begin their plunge, during which time detected signals may yield important information, 
as mentioned above. The lack of information about the dynamics of binaries before the plunge also hinders testing 
the validity of the initial matter and field configurations used in full GR calculations. While the assumptions defining 
the quasi-cquilibrium configuration are certainly reasonable, it is difficult to gauge how well such configurations agree 
with those evolved dynamically from larger separations. Also unknown in detail is the effect that errors in the initial 
configuration will have on the system as they propagate in non-linear fashion during the calculation. 

A further technical difficulty with full GR calculations results from computational limits, since numerical grids are 
currently constrained to have their boundaries lie within approximately half of a GW wavelength, within the near 
zone. This can induce possibly significant errors into the GW extraction process, which would ideally be performed 
by studying the behavior of the metric in the wave zone. 

A possible middle road, at least at present, is provided by the Conformally Flat (CF) approximation to GR, first 
described by Isenberg Q and developed in greater detail by Wilson and collaborators 0J ( n °te that their original 
version contained a mathematical error, pointed out by Flanagan |46j, accounting for spurious results with regard 
to "crushing" effects on NS prior to merger). This formalism includes much of the non-linearity inherent in GR, 
but yields a set of coupled, non-linear, elliptic field equations, which can be evolved stably. The first dynamical 
3-d calculations to make use of the CF framework were performed by Shibata, Baumgarte, and Shapiro |47], who 
created a PN variant by discarding some of the non-linear source terms in the field equations while retaining the vast 
majority of the non-linearity in the system. They found, among other things, that the maximum density of the NS 
is smaller for binary configurations than in isolation, and that NS in binaries have a higher maximum mass as well, 
strongly indicating that collapse to a BH prior to merger is essentially impossible. The first calculations to use the 
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full formalism were performed by Oechslin et al. |48|, using a Lagrangian SPH code with a multigrid field solver. 
Using corotating initial configurations, which are thought to be unphysical (since viscous effects are thought to be 
much too weak to tidally synchronize the NS fdH IHfT|L they found that mass loss during the merger is suppressed in 
relativistic gravitational schemes, as had been previously suggested based on calculations in PN gravity (FR3). 

While the CF formalism appears to be safer than full GR with regard to evolving the system stably, the grid-based 
approaches used so far suffer from many of the same problems faced by full GR calculations. Since the CF field 
equations are non- linear and lack compact support, approximate boundary conditions for the fields must be applied 
at the boundary of the grid, which can lead to errors in the solution. Additionally, very large grids are required to 
satisfy the tradeoff between large grid sizes, required to improve the accuracy of the boundary conditions, and small 
grid cell sizes, necessary to improve the resolution of the NS. Multi-grid techniques can be very valuable for such cases, 
but grid refinement techniques can also act as a source of numerical errors. Noting this, we have taken an entirely 
new approach to the dynamical evolution of binary NS systems with SPH, which requires no use of rectangular 3-d 
grids to solve for the metric fields. Instead, we use the LORENE numerical libraries, developed by the Meudon group, 
which are freely available online at http : //www . lorene . obspm . f r| (see [HlL hereafter BGM] and for a thorough 
description of the numerical techniques used in LORENE). These routines, which use spectral methods and iterative 
techniques to solve systems of coupled non-linear multidimensional PDEs, have been used widely to study quasi- 
equilibrium sequences of binary NS systems in Newtonian an d CF gravity (TG,[26|j5j|), as well as a number 

of other fluid configurations. The CF quasi-equilibrium solutions of Gourgoulhon et al. |26l hereafter GGTMB] and 
TG have been used as the initial configurations for the most recent full GR calculations of Shibata, Taniguchi, and 
Uryu, ;.43|, and will be used in this work as well. 

Our work here represents the first time that spherical coordinates and spectral methods have been used to study 
the dynamical evolution of binary NS systems in any gravitational formalism, including Newtonian gravity. It has 
long been known that these techniques are ideal for describing both binary systems with large separations as well as 
the merger remnants formed during the coalescence, since the spherical coordinates correspond much more naturally 
to the metric fields than rectangular coordinates can. Traditionally, rectangular grids have been used anyway, both 
because they are more widespread throughout computational fluid dynamics, but also because they are viewed as more 
robust. It has often been assumed that spherical coordinate field solvers may fail to calculate fields properly during 
the merger, when the matter can be described neither as two spheroids nor as one, but rather some combination of 
the two, with mass loss streams and other phenomena confusing the picture. We show here, however, that spectral 
methods can be used successfully in this regime, taking advantage of the multi-domain techniques of LORENE. This 
allows us to take advantage of the many advantages inherent to spectral methods: improved speed, vastly improved 
computer memory efficiency, and a coordinate system which allows for a natural treatment of the exact boundary 
conditions. 

The code we have developed to perform 3-d, relativistic, numerical hydrodynamic evolutions is well-suited for 
the study of a number of physical systems. While we focus here on merging binary NS, our code is capable of 
evolving essentially any binary or single-body relativistic fluid configuration. These include collapsing stellar cores 
and supermassive stars, and rapidly rotating fluid configurations. Modules currently exist to handle a number of 
physically-motivated EOS, including models with phase transitions and bosonic condensates, and we plan to extend 
our studies to include appropriate treatments of these conditions in the future. 

Our paper is structured as follows. In Sec.^J we summarize the theoretical basis for our new relativistic Lagrangian 
code, describing in turn the numerical methods used to implement both the dynamical equations of the CF formalism 
in SPH and the details of our spectral methods field solver. In Sec. IIIII we turn to the computational aspects of the 
code relevant to coalescing binary systems. We detail the choices we have made with regard to SPH discretization, the 
techniques used to convert between particle-based and spectrally decomposed descriptions of various quantities, and 
the coordinate transformations implemented to describe binary systems, merging systems, and the resulting remnants. 
In Sec. lIVI we report the results of several test calculations, including two well-known exact CF solutions as well as the 
evolution of quasi-equilibrium configurations in the absence of dissipative effects, producing circular orbits. In Sec.lVl 
we show our results from the calculation of a full NS binary coalescence started from the innermost quasi-equilibrium 
configuration, and followed through the merger and the formation of a merger remnant. We compare our results to 
previous work in the field, including full GR calculations using the same initial condition. Finally, in Sec. IVII we 
summarize our results and lay out some of the many classes of problems in relativistic hydrodynamics where our code 
may prove useful. 

II. NUMERICAL METHODS 

The CF approximation assumes that the spatial part of the GR metric is equal to the flat-space form, multiplied 
by a conformal factor which varies with space and time. Setting G = c = 1, as we will do throughout this paper, the 
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CF metric takes the form 

ds 2 = ~(N 2 - N.N^dt 2 - 2N l dtdx i + A 2 f tj dx l dx j , (2) 

where N is the lapse function, Ni is the shift vector, and A will be referred to here as the conformal factor (see Table|U 
for a comparison of our notation to those of [26l l45l I4H |48| , which are all based on the same exact assumptions) . The 
flat-space three-metric is denoted by /y. We follow the standard notation for relativistic tensors, denoting spatial 
three-vectors with Latin indices, and relativistic four-vectors with Greek indices. While the CF approximation cannot 
reproduce the exact GR solution for an arbitrary matter configuration, it is exact for spherically symmetric systems, 
and yields field solutions which agree with those calculated using full GR to within a few percent for many systems 
of interest 

As is typical for any relativistic calculation utilizing a polytropic EOS, there is a single physical scale which defines 
the units for the problem. We choose to define this scale for single-body test calculations (Sees. IIV Al and IIV B|) 
by dividing all mass, length, and time-based quantities by the gravitational (ADM) mass Mq of the object. For all 

binary calculations, we scale our results by the "chirp mass" of the system, A4 c h = p 3 ^ 5 Mf^ 5 , where M* is the total 
gravitational mass of the system at large separation, and p = M 1 M 2 /M t is the reduced (gravitational) mass. This 
quan tity is expected to be the most directly measurable physical parameter deduced from any GW observation (see, 
e.g., [53). For a binary system consisting of two Mq = 1.4 Mq NS, the chirp mass is M. c h — 1-22 Mq, which yields 
characteristic distance and time scales R = 1.80 km and T — 6.00 x 10~ 6 sec, respectively. To compare our results 
with the equivalent relativistic model (run M1414) of 0], who use the total gravitational mass of the binary at large 
separation as their unit, one can divide our quoted masses and radii by a factor of 2 1 ' 2 = 2.30. To convert our time 
evolution results into the time units used in [43| , which are defined in terms of the initial binary orbital period, divide 
our time units by a factor of 443 (their unit Pt=o = 2.66 ms for two NS each with Mq = 1.4 Mq). 



A. CF Smoothed Particle Hydrodynamics 



In what follows, we will assume that the stress-energy tensor of the fluid is that of a perfect fluid, with 

Tpu = (p + pe + P)u u ii„ + Pg^, (3) 

where p is the rest-mass density, e is the specific internal energy, P the fluid pressure, and u n the four-velocity of the 
fluid. For our initial data, we assume a polytropic EOS P = kp r , with constant values for k and the adiabatic index 
r, and throughout the calculation we assume that the evolution is adiabatic, such that P = (T — l)pe. The maximum 
infall velocity of the two stars relative to each other during our dynamical calculations was found to be o m = 0.06c, 
whereas the sound speed at the center of the NS is initially c s = 0.85c. Since the sound speed for our EOS depends on 
density such that c s cx p 1 ^ 2 , we expect that the only supersonic fluid flows will occur in the very tenuous outer regions 
of the NS. During the merger, we do expect that low-density matter from the inner edge of each NS will cross over to 
the binary companion at high relative velocity (Sv ~ 0.3c), but the velocity field is dominated by the circular motion, 
rather than a converging flow. All motion within the cores of the NS should remain strongly subsonic throughout the 
merger process. 

Using our CF metric and stress-energy tensor, Eqs. [21 and [21 the Lagrangian continuity equation is given by 

^+^ v y = o, (4) 

dt 

where the rest mass density is defined by 

p* = Nu°A 3 p = ln A 3 p, (5) 

the Lorentz factor of the fluid is defined as 

In = Nu°, (6) 



and the physical velocity is given by 



uW = N* + -^. (7) 



Note that u° is the timelike element of the covariant 4- velocity; all other numerical superscripts refer to exponents. 
Following the SPH prescription, this conservative form of the continuity equation allows us to define a set of particles, 
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each of which has a fixed mass m a and a well defined velocity given by dxi/dt = Uj. For each particle "a", we 
also define a "smoothing length", h a , which represents the physical size of the particle. SPH particles do not have 
delta-function density profiles; rather, each particle represents a spherically symmetric density distribution centered 
at the particle's position x a with compact support. This density distribution, determined by our chosen form of 
the smoothing kernel function, is second-order differentiable, and drops to zero at a radius equal to two smoothing 
lengths from the particle's position. For each particle, we define a set of neighboring particles by the condition that 
all particles whose centers fall within a given particle's compactly supported density distribution are its neighbors. To 
determine the proper smoothing length for each particle, we define an ideal number of neighbors for each particle, iV/v, 
and we use a relaxation technique to adjust h a after every timestep (as described in detail in, e.g., FR1 and HU). We 
note in passing that "neighborhood" is not a reflexive property; we handle this through the use of a "gather-scatter" 
algorithm (see [59j for details). 

The primary advantage of the SPH method over traditional grid-based codes is that fluid advection is handled in 
a natural way, such that one can define the edge of a fluid distribution without recourse to artificial "atmospheres" 
or other tricks necessary to prevent matter from bleeding into the vacuum. Particles simply follow their trajectories 
in the fluid flow. As such, SPH is extremely computationally efficient, since all numerical resources are focused 
automatically on those regions containing matter. Because of this, SPH also allows for high spatial resolution. The 
primary disadvantage of SPH compared to shock-capturing grid-based codes is in the lower resolution of shock fronts, 
which, as we discuss below, is unimportant here. 

All hydrodynamic quan tities can be defined using standard SPH summation techniques over each particle's neighbor 
list (see, e.g., |58l l60lT6lT | for a thorough discussion), with the rest mass density p* taking the place of the standard 
Newtonian density. Thus the rest-mass density for particle "a" can be defined via SPH summation over a set of 
neighboring particles denoted by "b" as 



(P*)a = rn b W ab 



(8) 



where W a b is the Wa smoothing kernel function for a pair of particles first introduced by Monaghan and Lattanzio 
|62| , used in FR, |58| and many other implementations. The momentum equation is given by 



(9) 



where the specific momentum is defined by 



= huj 



(10) 



and the specific enthalpy is defined as 



h = (1 + e + P/p) = 1 + Te. 



(11) 



In this expression and throughout this paper, covariant derivatives are associated with the flat-space metric. In the 
absence of non-adiabatic artificial viscosity terms, the energy equation merely implies that the value of k in the EOS 
remains constant. In our calculation, we use p* and Ui as the basic hydrodynamic variables (in addition to our 
uniform value of fc). To find u , which enters into the momentum equation, we take the normalization condition for 



the 4-velocity, 



= — 1, and find 
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(Nu 



0\2 



UiUi UjUj 
1 + -7T- = 1 + 



A 2 



A 2 



Tkp 



r-i i 



(7n^) 



r-i 



(12) 



which can be solved implicitly in terms of the density and the field values. 

To solve the field equations of GR, we need to fix the slicing condition for timelike hypersurfaces. Following the 
standard approach to the CF formalism, we find that the extrinsic curvature tensor is given in terms of the shift 
vector as 



JC 3 = 



1 



2NA 2 



V 3 N* - -f j V k N k 
3 



(13) 



when we assume the maximal slicing condition (TrK = 0). To lower the indices of the extrinsic curvature tensor 
and other spatially-defined quantities, we use the conformal flat-space metric (i.e., Kij = A^Sn-SjiK^ in Cartesian 



6 



coordinates, where 8ij is the Cartesian flat-space metric). Combining the maximal slicing assumption with the 
Hamiltonian constraint yields a pair of elliptic equations for v = \n(N) and (3 = ln(A/V) (GGTMB), in the form 



VfcV'v = 4nGA 2 {E + S) + A 2 K i3 K l: > -VtisV^, (14) 



V fe V fe /3 = 4irGA 2 S + -A 2 K l2 K^ - -(V^VV + V,/3V l /3), (15) 



where the matter energy density and the trace of the stress energy tensor are given by 

E = ln 2 {ph)-P, (16) 

-Y 2 - 1 

S = ,n . {E + P) + 3P. (17) 

In 

Note that Eqs.^]andEl are algebraically equivalent to the field equations found in other papers on the CF formalism, 
although most groups have typically solved the corresponding Poisson-type equations for ip = \f~A and Nip. Finally, 
the momentum constraint gives us the equation for the shift vector, 

V^'VjiV* + iv'Vj-iV 3 ' = -16ttG— (E + P)ui + 2NA 2 K ij V j (3f3 - Au). (18) 

3 h~f n 

Since the CF formalism is time-symmetric, and we can define several conserved quantities. The total baryonic mass, 

M b = / p*d 3 x, (19) 



is automatically conserved in our SPH scheme, since we use the rest mass density to define particle masses. The total 
angular momentum of the system can be defined as 



Ji = £ijk J p*x 3 u k d x. (20) 

Finally, the ADM mass is given by 

M a = J PADMdPx; Padm = A 5 / 2 \E + -L^K^rA , (21) 

It is important for numerical reasons to note that the two terms that make up the ADM mass density padm have 
different behaviors: the contribution from the matter energy density, 



Pi 



A 5/2 E, (22) 



has compact support, and is non-zero only in the presence of matter, whereas the term involving the extrinsic 
curvature, 

p^^A^K,, (23) 
extends throughout all space with power-law fall-off at large radii. 



B. The Spectral Methods Field Solver 

The spectral methods techniques we use to solve the CF field equations, Eas. 1141 1151 and fTTH discussed in great 
detail in |52j, provide a number of important advantages not present in traditional grid-based approaches. First and 
foremost is their speed and computational efficiency. Finite differencing schemes typically require 3-d grid sizes of 

10 6 elements 0, 48]. In contrast, GGTMB show that spectral method field solvers can be used to construct field 
solutions yielding ADM masses and angular momenta convergent to within 10~ 4 and satisfying the virial theorem 
to the same level using only 3 grids of size 17 x 13 x 12. These grids, extremely small compared to those used in 
Cartesian multi-grid solvers, result in a great increase in speed. Additionally, the use of spherical coordinates allows 
for a more natural treatment of boundary conditions. In the approach we use, taken from GGTMB and summarized 
here, space is described using spherical coordinates, split into a set of nested "computational domains" . The outermost 
domain can be compactified by rewriting the field equations in terms of 1/r, allowing us to impose the exact boundary 
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conditions at infinity, rather than "fall-off" boundary conditions which approximate the behavior of the fields at the 
edges of a rectangular grid. As such, we avoid the classic tradeoff between a grid with large grid spacing, which yields 
accurate boundary conditions but poor spatial resolution of the matter source, and a grid with small spacing, and the 
opposite concerns. 

Combining the LORENE methods with particle-based SPH requires some small but significant changes from the 
previous approach described in detail in Sec. IVA of GGTMB. As in that work, we construct a set of three compu- 
tational domains around each star to evaluate the field equations. The innermost domain has a spheroidal topology, 
with a boundary roughly corresponding to the SPH particle configuration's surface, as described in Sec. IHI Al The 
other two domains cover successively larger regions in radii, with the outermost domain extending out to spatial 
infinity, as shown in Fig. The field equations are solved in each domain and the global solution is obtained by 
matching the function and its first radial derivative at each boundary. Appropriate boundary conditions at radial 
infinity are also imposed. All fields can be described in one of two complementary representations, either in terms 
of their coefficients in the spectral decomposition, or by means of their values at a set of "collocation points" . The 
coordinate representations for these points are defined such that the origin of the system describing each NS is located 
at the position of the star's maximum density, The collocation points are spaced equally in both sin# and <f>, as 
required by the angular component of the spectral decomposition. The radii of the points are determined from the 
collocation points of the Chebyshev polynomial expansion used for the radial coordinate, as described in BGM and 
GGTMB. For all calculations described in this paper, these domains consisted of a 17 x 13 x 12 grid, in terms of radial, 
latitudinal, and longitudinal directions, respectively, an acceptable trade-off between speed and accuracy, based on 
the description above. 

A key feature of the LORENE libraries is their handling of binary systems in a straightforward and natural way. 
This involves "splitting" the source terms of the field equations into two distinct components, each of which is centered 
on one of the stars in the binary. Since the field equations in this case are non-linear the split cannot be performed 
uniquely; the fields present around one star cannot be determined from its source terms only. Rather, this method 
seeks to minimize the contribution of one star to the fields around the other. In practice, both field variables and 
hydrodynamic source terms can be broken down, as shown in Eqs. 80-88 of GGTMB, such that 

v = v<\> + v<2> = v<i> + v <2 ^i> = t / <i^2> + ^<2>, etc. (24) 

where quantities labeled < 1 > and < 2 — > 1 > are defined at the collocation points of star 1, and < 2 > and 
< 1 — > 2 > at those of star 2. The autopotentials of each star, v < \ > and i^<2>, are primarily generated by matter 
from the star itself, while the "comp-potentials" , ^<2^i> and ^<i^2> are primarily generated by the other star. It 
is this conversion of fields between the two sets of coordinates which represents the greatest amount of numerical 
effort during a calculation. In practice, we attempt to minimize the magnitude of the comp-potentials since they are 
centered around the other star and not as well described by spherical coordinates. The detailed description of how 
these quantities are defined and calculated can be found in Sec. IVC of GGTMB. 

III. NUMERICAL TECHNIQUES FOR COALESCING BINARIES 

Integrating the LORENE library routines into an SPH-based Lagrangian code introduces a number of rather subtle 
numerical issues. The simplest of these is deciding the shape of the innermost computational domain's surface. The 
closer the boundary of the innermost domain lies to the surface of the fluid, the more Gibbs effect errors are minimized, 
so long as the surface is sufficiently smooth and convex everywhere. Gibbs effects, which are common to all spectral 
decomposition techniques, result from the attempt to describe a non-differentiable density distribution as a weighted 
sum of smooth functions; they are relevant here when we attempt a spectral decomposition of a region of space 
where the density drops to zero within the boundary (see BGM for more details). Unfortunately, the smoothness and 
convexity conditions do not necessarily apply to the region in space where the SPH density is non-zero, since a single 
SPH particle being shed from the star can greatly affect the non-zero density boundary, even though it may represent 
an insignificant amount of mass. As a result, we have implemented an algorithm that attempts to take the middle 
ground, defining a boundary for the innermost domain which encloses as much of the mass as possible, in such a way 
that smoothness and convexity are guaranteed. The second domain is bounded by a sphere of radius twice that of 
the outermost point of the first domain, and a third domain extends to spatial infinity. 

Once the configuration of the computational domains is determined, there are several choices which need to be 
made with regard to the most accurate and efficient way to calculate various terms in the evolution equations. Some 
hydrodynamic quantities, such as the rest mass density p*, need to be defined for each of the SPH particles. Here, 
we use N ~ 100, 000 SPH particles for each run, which was found to be sufficient for achieving numerical convergence 
of the GW signal to the ~ 1 — 2% level in our studies with post-Newtonian SPH 35]. Other quantities, such as 
those appearing in the source terms for the field equations, need to be defined at every point among the 17 x 13 x 12 
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spheroidal grids of collocation points for each of the three domains. Compared to solutions computed using larger 
grid sizes, we find that these agree with significantly larger grids to within ~ 0.1% for the value of the shift vector, 
and even better for the values of the lapse function and conformal factor. 

In general, however, many quantities do not need to be defined both ways, so long as a full set of thermodynamic 
variables is known in both representations. Details about which quantities are used in each representation are given 
below. Briefly, it is most efficient to perform the majority of our algebraic operations on quantities defined at 
collocation points, reading in and exporting back as small a set of parameters as possible to the full set of SPH 
particles. While reading quantities from SPH particles to collocation points is relatively quick, requiring an SPH 
summation at every collocation point position, the reverse process is much more involved. To calculate the value of 
a quantity known in the spectral representation at SPH particle positions requires performing a sum over all spectral 
coefficients with the weights appropriate to each particle position, consuming a great deal of time. 



A. Binary systems 

To construct the initial SPH particle configuration for a binary NS evolution, we use the quasi-equilibrium irrota- 
tional models of TG, which are publicly available at http://www.lorene.obspm.fr/data/ These models describe 
the complete 3-dimensional structure of both the field values and hydrodynamic quantities at every point in space, 
in terms of a spectral decomposition. Specifically, we take the results of their irrotational run for stars of equal mass 
and equal compactness GM/Rc 2 = 0.14 with a T = 2 polytropic EOS. Each star has a baryonic mass in isolation 
given by GM B /Rc 2 = 0.1461. 

To convert these models, which are stored in the coefficient basis, into spatially defined particle-based quantities, 
we first lay down a grid of SPH particles in a hexagonal close-packed (HCP) lattice with constant lattice spacing and 
particle smoothing length. This grid is then treated as if we reflected around the z — plane, to take advantage of 
the vertical symmetry inherent in the problem. Each particle is treated as if it were really two particles for all SPH 
summations, one located above the z = plane, the other an equal distance below, both with half the true mass of 
the "real" particle. Since the vertical symmetry is enforced on a particle-based level, we solve all our field equations 
only for vertical angles < 8 < tt/2, and reflect the solutions for all points below the plane. The mass of each particle 
is initially set to be proportional to the density at the particle's position according to the quasi-equilibrium model. 

Next, using an iterative process, we calculate the SPH expression for the density of each particle, using Eq. |H1 and 
adjust each particle's mass so that the SPH density matches the proper value from the initial model, stopping once 
the maximum difference for any particle is less than 0.25% of the star's central density. Particle velocities are assigned 
to match the quasi-equilibrium model's velocity field, as are all other thermodynamic variables. Finally, we advance 
the velocities by half a time-step, using the same methods as in a standard iteration loop (see below), since we use a 
second-order accurate leapfrog algorithm (described in detail in HU). 

During each iteration, the first step is always the calculation of each particle's neighbor list, and the associated 
SPH forms for the density and other hydrodynamic terms. Once this is done, we perform a Euclidean transformation 
on our coordinates into a new frame (denoted by primed quantities), whose origin is defined to be the system center- 
of-mass, with the center-of-mass of each star lying on the x-axis, making sure to transform the positions, velocities, 
and accelerations for each particle. In terms of the inertial frame coordinates of the NS centers-of-mass, x\ and x%, 
the transformation is given by 

^ = tan- 1 l^i, (25) 
x 2 - x\ 

_ x! + x 2 _ yi + 2/2 , 0RA 

xcoM = 2 > VCoM = — 2 — ' 

x'(x,y) = (x - x C om) cos 4>' + {y - ycou) sin 4>\ (27) 

y'(x, y) = (y- ycoAi) cos cf>' - (x - x Co m) sin 0', (28) 



R 



'1,2 



(±-,0, 0);R= y/(x2-xi) 2 + (y2~yi) 2 , (29) 



This transformation is shown in the upper left panel of Fig. [21 The next computational task is defining the shape 
of the innermost computational domains, calculated for sets of rays equally spaced in and <p, as measured from the 
center-of-mass of each star. We denote these surface functions r(9 q , <fi q ), where the "q" subscript, taking the value 1 
or 2, refers to angles measured in the primed frame of Eas. I25H29I outward from the center-of-mass of star q. These 
surfaces are used to determine the position of the collocation points, which in turn are used as the basis for the 
entire spectral decomposition. While it is easy to find the point along any ray where the SPH density drops to zero, 
rsPH(0 q ,4> q ), we have found that such a set of points leads to unacceptable results from the field solver, especially 
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with regard to the convexity of the matter distribution. If we ignore the density contributions of all particles whose 
density falls below some fixed value, say p m i n = 0.0001, the resulting surface function fspH{0 q ,4> q ) is typically much 
more regular, but still is not an ideal choice, since we still cannot guarantee convexity. It should be noted that the 
field solver is entirely capable of handling matter sources which lie outside the innermost domain, although it does 
work best in the case where the surface of the matter matches the domain boundaries closely. For this reason, we 
restrict the shapes of the innermost domains to triaxial ellipsoidal configurations, oriented along the principal axes of 
the moment of inertia tensor for each star. The growing misalignment between the stars and the (rotated) coordinate 
system is a well-known effect, reflecting the tidal lag that develops in close binaries as matter tries to reconfigure itself 
in response to the rapidly changing gravitational field /citeLRS5. Thus, for each star, we calculate the angle Q q such 
that 



= i tan" 



21 



xy 



L U!) 



where iy = ^ a m a (x' a — x' q )i(x' a — x' q )j, and define our surface function r(Q qi (j) q ) such that 

;21 -1/2 



r(0 q ,(p q ) 



l> 2 



where 



(30) 



(31) 



'q, Vq 



) = 



z{6q,<Pq) 



sin 6 q (cos $ q cos i 
sin # g (cos$ 9 sin^ 
cos 9„, 



- sin <fr q sin cj> q ) 
sin <i> g cos cj> q ) 



(32) 
(33) 
(34) 



and a, 6, and c are the axis lengths of the ellipse. 

We have found it best to fix the axis ratios of the ellipse by computing the maximum extent of the surface defined 
by r(6 q , 4> q ) in each principal direction, such that 



a = max | t S ph( 
b = max | r S pH( 
c = max | r S PH{ 



'q, i 



%(8q,<f>q)\, 
V{0q,<t>q)\, 
Z{0q,(t>q)\- 



(35) 

(36) 
(37) 



Since this prescription can lead to some particles with p a > p m i n falling outside of the surface, we multiply the 
distance in all directions by the smallest factor Fq required to encompass all SPH particles whose density is greater 
than p m in, typically leading to an increase in linear size of no more than 2%, such that a = ao-Fbj b = boFo, c = cqFq, 
and 



F 




(38) 



A typical surface fit is shown in the upper right panel of Fig. [5] The outer two domains are defined in terms of these 
new coordinate systems as well, to allow us to match field values and their derivatives at the boundaries. 

Using the collocation points derived from these surfaces, our next task is to calculate the value of the field equation 
source terms at these points. Since calculating E and S for each SPH particle from Eas. 1161 and 1171 would require 
a great deal of ultimately needless algebraic work, we do not calculate the SPH expressions for these quantities at 
collocation points. Instead, at every collocation point, we calculate the SPH expression for the rest mass density, the 
"rest pressure" P* = kp^, and the density- weighted average velocity Uj. Using these values, as well as the field values 
from the previous iteration, we calculate 7„ from Eg. 1121 and then E and S at every collocation point, and proceed 
to solve the field equations, using the iterative techniques described in GGTMB. Typically, we require « 50 iteration 
loops to achieve a solution for which no field value varies by more than 1 part in 10 9 from one iteration to the next. 

Once we have solved the field equations, we use the spectral expansion of the fields to calculate as much as we can 
of the terms in the force equation, before reading off the values at every particle position, which takes a significant 
amount of time. In practice, we use the spectral decomposition to calculate the prefactor for the pressure force term, 
NA 3 , the vector sum of the terms involving derivatives of the conformal factor and lapse function, the nine first 
derivatives of the shift vector, and the radiation reaction terms. While it would be faster to calculate the force term 
involving the derivative of the shift vector completely in the spectral basis, we have found it to be inadvisable. This 
term alone is linear in the velocity, and it is inconsistent to use an averaged velocity in this term on the RHS of 
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the force equation to calculate the rate of change of each particle's individual velocity on the LHS. Thus, denoting 
terms calculated within the spectral basis and exported to particle positions by "sb" , and those calculated using SPH 
techniques only by "SPH" , Eq. is truly evaluated as 



dui 
~dt 



= [-NA 3 ]. 





+ 


p* 


SPH 



-ViA - hj n S7 l N 



~ [W l N : >} sb [u 3 }sPH. 



(39) 



After calculating the forces on each particle and advancing the velocities by a full timestep, we are still left with the 
task of recomputing Eq. [T] which relates v l and Uj. Since the sources for the field equations are velocity dependent, 
and we have just advanced u, we rerun the field solver with the new values of Uj, and compute Eq. \7\m the form 



v l = [N% b + 



1 



A 2 hu n 



[Ui\sPH- 



(40) 



Finally, we record the GW strains, ADM mass and system angular momentum, after rotating all positions, velocities, 
and accelerations back to the inertial frame by means of the inverse Euclidean transformation to the one at the 
beginning of the iteration 

There are several different ways to calculate the ADM mass numerically, all of which should be equivalent to Eq. 
yielding an important check on the code. First, we calculate the system's ADM mass using by taking a surface integral 
at spatial infinity (see Eq. 65 of GGTMB), 



Mo = ~i d l A 1 ' 2 dS l 

^ Joo 



(41) 



This quantity can be compared to the particle-based expression, with two important caveats. Since the extrinsic 
curvature contribution to the ADM mass, Eq. 1231 does not have compact support, there is no way to convert the 
integral into a sum over particles that have a different spatial extent. Second, noting our concerns about exporting 
large numbers of terms from the spectral representation to particle positions, we perform some of the algebra involved 
in determining the ADM mass in the spectral representation. In the end, the particle-based expression for the ADM 
mass becomes 



M = (J p 2 d 3 x) sb + y^m a 



(42) 



and we define p = pi/ p* as the ratio of ADM mass density to rest mass density, calculated using LORENE techniques 
at all SPH particle positions. 

In a similar fashion, there are two ways to calculate the system angular momentum. We check the behavior of J in 
the spectral basis (see Eq. 67 of GGTMB), 



167T Jo,, 

as well as the SPH expression for the angular momentum, 



m a x J Uk. 



(43) 



(44) 



Since the CF formalism is time-symmetric, the dissipative effects of gravitational radiation backreaction have to 
be added in by hand, just as they are for PN calculations. Previous PN calculations of binary NS systems (FR and 
[53) have typically employed the exact 2.5PN formulae introduced by (Zjj to describe lowest-order GW losses from 
the system. Unfortunately, those equations are not applicable to CF calculations, since they are written in terms of 
fields defined in the P N ap proximation that differ from those defined in the CF approximation. For this work, we 
follow the approach of |45|. using the slow- motion approximation to estimate the radiation reaction potential of the 
system. While this method contains some obvious flaws, most obviously the fixed spatial dependence of the radiation 
reaction potential, it does yield a backreaction force which is quantitatively correct in overall magnitude. These 
approximations should not affect our calculations to a large degree. While the infall velocity of the binary prior to 
plunge is driven by the GW backreaction, a different regime occurs after dynamical instability sets in. During this 
period, the evolution is almost completely hydrodynamic in nature [S8| . While the chosen GW backreaction treatment 
may affect the final mass and angular momentum of the resulting merger remnant, it will play only a secondary role in 
the detailed evolution of the fluid configuration, since GW backreaction become less important during the coalescence. 
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Prom Eq. 51 of |4q|. the radiation reaction force in the slow- motion approximation is given by adding a term, 

a t:r eac = N 2 h U °V lX , (45) 

to the RHS of Eq. 03 We define the radiation reaction potential x m a similar but slightly different way than their 
Eq. 56, such that 



X 



^x k x Q\J, 



similar to the approach taken in |48j |. As they do, we define the quadrupole moment as 



I hi 



STF 



PADMXkXld X 



(46) 



(47) 



noting that our "padm" corresponds to many other authors' S^,, or some multiple thereof. The expression "STF" 
refers to the symmetric, trace-free component of the tensor, which is linked to the gravitational radiation production 
in the quadrupole limit. As noted previously in our calculations of the system's ADM mass, we can evaluate the 
contribution from p\ using standard SPH summation techniques, but the second term can only be found using 
LORENE integration techniques. Thus, 



Qki = Qi + Q 2 = STF 



^m a p a x k a x l a ) + (^J d 3 xp 2 x k x l 



sb 



(48) 



where Q\ and Q2 reflect the contributions from p\ and p 2 , respectively. Using SPH, we can take the first time 
derivative of the former half, so long as we ignore the Lagrangian derivative dpi/dt, which should be essentially 
negligible during our calculations. We find 



(Q 



X)kl 



STF 



Pi{x k vi + xiv k )d 3 x 



(49) 



To calculate the rate of change of the extrinsic curvature, we assume that the time variation in the tensor is due solely 
to the orbital motion (rather than an overall change in magnitude of the tensor components in a corotating frame), 
which yields 



— {Qilyy 



2uj(Q 2 ) X y, 
{Q2)xx 



2uj 



(Qz)yy 



(50) 
(51) 
(52) 



where the factor 2u> reflects the fact that the quadrupole tensor makes two cycles during every orbital period. To 
calculate the fifth time derivative of the quadrupole tensor, we use the same technique with both components of the 
tensor, finding 



[5] 



16w 4 Qa 



(53) 



where in all cases the system's instantaneous angular velocity is calculated as the ratio of the angular momentum to 
the moment of inertia, 



= J2a m a (x a (v y ) a - y a (v x )g) 



(54) 



which holds exactly in the quadrupole limit for synchronized binaries. 

Calculating the GW signal and energy spectrum is a more straightforward task which can be done after the 
calculation is finished. We calculate the GW strain in two independent polarizations for an observer located at a 
distance d from the system perpendicular to the orbital plane from the lowest-order quadrupole expressions, 



dh-i 
dh y 



[2] - 



[2] 

yy 



2Q 



■nr 



(55) 
(56) 
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where we calculate the second time derivatives of the quadrupole tensor by numerically differentiating the results 
from Eq. 1491 In terms of the Fourier transform of the quadrupole moment, 



= 27ri/t n [2] 



Q[i(t)dt, 



(57) 



where few = 2forb is the G W frequency, the GW energy spectrum is computed as 
dE 



dfew 



(Q 



2] 

XX 



(Q [ * ] y) 2 + mr+m) 



)[2]\2 



(58) 



B. Merging systems 

As the stars in the binary system spiral inward, they reach a point where the density distributions begin to overlap 
as matter from the inner edge of each star falls onto the surface of the other. Our field solver can handle this situation, 
since it does not assume that the matter sources are spatially distinct, but the surfaces required to envelop the particles 
from each star would become poorer and poorer fits to the two NS. Noting this, we alter the approach described above 
in a number of ways when the particles first cross through the inner Lagrangian point at the center of the system, in 
such a way that the surfaces we define for each star always remain smooth and still do an acceptable job of describing 
the true density distribution in a meaningful way. 

Once the binary separation shrinks sufficiently, matter streams from the inner edge of each NS toward the other 
NS, flowing along the surface of the companion. These counter-streaming, low-density flows lead to the formation 
of a vortex sheet (FR3). Mass transfer typically occurs before the triaxial ellipsoidal surfaces used to define the two 
stars overlap, since particles crossing from one star to another generally fall beneath the density cut used to define 
each surface. To account properly for the star to which each particle is bound, we declare particle "a" a member of 
the first star if x' a < and a member of the second star if x' a > 0, at least for the purposes of defining each star's 
center of mass and tidal lag angle, as described in Eqs. I25H38I 

This approach entails further alterations once the ellipsoidal computational domains from each star begin to overlap. 
Using the j/'-axis as the dividing line between the two stars is fine for determining the center-of-mass, tidal lag angle, 
and ellipsoidal surface for each star, but it is inappropriate to draw a fixed line at x' = when calculating field 
equation source terms. These would induce a sudden density drop from finite density at small negative values of x' 
to zero density at positive x' values (for the first star, vice versa for the second), leading to large Gibbs effects. It 
is equally inappropriate to count one particle as a member of both stars, since we would end up double-counting its 
density contribution to the field equations. Instead, we introduce a weight function f a = f(x' a ), for each particle in 
the overlapping region in such a way that / = at the surface of star 1, / = 1 on the surface of star 2, and < / < 1 
in a spatially differentiable way within the overlap region. To do so, we define /i and /2, the fractional squared 
distance outward a point lies from the center of each star to the surface. For the first star, 

(X') 2 (Y') 2 (z') 2 
/,(x') = i-|L + i^L+l^_ , (59) 
a z b z c l 

where 

X' q = x q cos§ q + y'sin$ 9 , (60) 
Y' q ee -x g sm$ q + y'cos$,j, (61) 

are the rotated coordinates used to define the (tidally lagging) ellipsoidal surface of the star, and x q (x') = x'~x' q is 
the distance in the x'-direction from the center-of-mass of star q. A picture of these quantities is shown in the bottom 
left panel of Fig. [21 In terms of f q , we define our overlap function f a for each particle such that 

fa = f@ a ) = {1 _ f [ )3 / 1 ( [_ h)3 - ( 62 ) 

For the first star, source terms in the field equations are evaluated as 



G°*)i = ^2m a f a W a , 

a 



(63) 
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( p *)i = r^%W fc [(p*)i + (p*) 2 ] r , (64) 

[p*)l + (P*)2 

K)l = ~ ^ JTTr > (65) 

La m aJaW a 

where W a is the smoothing kernel function evaluated between the collocation point and overlapping particles. Com- 
plementary expressions hold for the second star by substituting (1 — f a ) for f a . 

This approach allows us to calculate the fields properly using spectral methods well into the merger, but a final 
modification is necessary to bring us to the point where the density profile of the matter can better be described by 
viewing it as a single object. When the inner edge of one star overlaps the center of the other star, the density profile 
typically becomes bimodal, leading to spurious results from the spectral expansion. To guard against this happening, 
we use a relatively simple approach. If the surface of one ellipse extends more than halfway from x' — toward the 
center-of-mass of the other star at x' q — ±R/2, we linearly rescale all surface points that lie across the y'-axis. Thus, 
defining x q (9 q , <j) q ) = r(8 q , 4> q ) sin0 o cos^> in accordance with our previous notation, if the surface of either star extends 
to a maximum value x q[max = max(|cc |) > 3i?/4, (or in other words, if max(|r(# , <p q ) sin0 o cos0 9 — x' q \) > -R/4,) we 
adjust all points on the other side of the y'-axis, yielding 



(a ^J,^ x q:old (0 q ,<j> q )-R/2 



for all points with x q (9 q , 4> q ) > R/2. (66) 



^q-.raax \ w qi 1 

To evaluate the weight function, we adapt the x'-dependence correspondingly, such that for particles with x q > R/2 

fa:new(x q ^) — fa:old{^ ' ^<y); (67) 

where the rescaling factor k = R/2 + (x q — R/2) Xg ' m ^ i R ^ 2 . The bottom right panel of Fig. [3 demonstrates this last 
coordinate transformation. 

Eventually, the system will reach a point where it can no longer be properly described as a binary, and our field 
solver fails to converge. Before this happens, we reach a point where we can describe the system as a rapidly rotating 
single star, using all the methods described above for computing the evolution, but now assuming that all particles 
comprise the same star. We have found that our results are independent of the exact moment at which we make 
this conversion. In the following discussion, we will show the results for runs where we perform the conversion at the 
earliest possible time for which the field solver will converge to a solution when the matter configuration is treated as 
a single star. 



IV. TEST CALCULATIONS 



We have performed several tests to check the accuracy and numerical stability of our code, for both single-star and 
binary systems. We have studied the behavior of the code for spherically symmetric problems for which the exact 
field solution can be calculated semi-analytically: the Oppenheimer-Volkov solution for a static spherical star and the 
collapse of a pressureless dust cloud initially at rest. We have also calculated the dynamical evolution of the binary 
quasi-equilibrium models calculated in TG, without the inclusion of radiation reaction forces. A summary of all our 
calculations, including those used for testing the code, can be found in Table ITT1 



A. The Oppenheimer-Volkov (OV) Test 

Since the CF formalism is exact for spherically symmetric systems, which can always be described in isotropic 
coordinates, it is fair to expect that any working code should be able to reproduce the well-known Oppenheimer- 
Volkov solution exactly, noting that the traditional form of the OV solution needs to be rewritten into CF coordinates. 
With p' = p(l + e) as the total energy density (including both the rest and internal energy densities), the OV equations 
are typically written 

dm 9 , .„„. 

— = 47TFV, 68) 
dr 

dP (pf + P)(m + AnPf 3 ) 

dr r z — Zmr 

d$ (m + inPf 3 ) 
dr f 2 — 2mf ' 
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for an interior metric in the form 

ds * = _ e ^ dt 2 + ( ! _ ^1 ) df 2 + f 2 dn2) 



2m 



and exterior metric in the Schwarzchild form 



where the star's total gravitational (ADM) mass Mq = m(f s ), and f s is the Schwarzchild coordinate radius of the 
stellar surface. A quick comparison with the CF metric, Eq. [3 shows that N = e* inside the star, but the conversion 
between A and $ requires more care, since we have to determine the relationship between the two coordinate radii 
r(f). The simplest way to determine the change of coordinates is to solve for the values of r s and f s , the radii of the 
stellar surface, using the asymptotic behavior of the exterior solution (following the same logic used in exercise 31.7 
of H3|)- Comparing the radial and angular parts of the metric, we find 

A 2 dr 2 = (l-^.) 1 dr \ (73) 



A 2 r 2 = f 2 . (74) 

Dividing and taking a square root, we find 

dr df 



(75) 



Integrating yields 



and we find 



kr 



r y/f(f - 2M Q ) 

lnr + k = 2 In (Vf + y/f - 2M ) , (76) 
(y/f + \/f - 2M ) 2 = 2f - 2Af + 2y/ f(f - 2M ). (77) 



The asymptotic behavior at infinity indicates that we must have k = 4, so our final expression for r(r) takes the form 

r = ~ (f - M + y/f(f - 2M )) . (78) 
For a star with Schwarzchild surface radius f s , we find that the CF surface radius is given by 

r s = \{r s - M + y/f 3 (f s -2M ) S j , (79) 

the surface value for the conformal factor is 

A s =- = 2f s (80) 

r f s - M + y/f(f-2M Q ) 

and the lapse function at the surface is given by 



(SD 

V 

To solve for the interior metric, we add an equation to the OV set to take into account the different radial schemes, 
defining a scale free conformal radius satisfying a boundary condition limf—,0 Tq = f whose radial behavior is given by 

^ = , r ° ■ (82) 
dr y/f(f - 2m) 
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This yields an expression for ro(r) which is defined up to an arbitrary multiplicative constant fe, such that tq — k-f(f), 
where / is determined from the mass distribution. To find r(f), we merely set k — r s j (r ) s , which implies r = fr s /(r ) s 
and determine A from Eq. 1741 

To test the code, we ran three calculations whose results could be compared with well-known semi-analytic solutions. 
First, we constructed an equilibrium model for an isolated NS, whose density profile was given by an OV solution 
for a r = 2 EOS, with unit ADM mass, and a conformal radius r s — 6.874. The solution has a total baryonic mass 
Mb = 1.066, and an areal radius f s — 7.913. We started run OV1 by taking this solution as an initial condition, 
and used it to test the overall stability of the code for equilibrium configurations. Next, we constructed a similar 
model with the same mass and EOS, but scaled to an initial radius 10% larger. Run OV2 is a dynamical calculation 
started from this initial condition using our standard evolution code, enabling us to study the oscillations around 
our equilibrium solution. Last, we took the non-equilibrium configuration from run OV2, but added a "relaxation" 
drag term of the form —Ui/t re i ax to the RHS of the force equation, Eq. [51 with t re i ax /Mo — 7.9. This provides an 
overdamped force for run OV3 since the dynamical timescale tp/Mo = (Gp)~ 5 = 18. All three runs were followed 
until t/Mo = 120, corresponding to 6.7 dynamical times for the NS, and in all cases radiation reaction was turned off. 

In Fig. we show the radial profiles of the lapse function and conformal factor at t/Mo = 100 for runs A and C, 
along with the correct OV solution. It is no surprise that run A has essentially remained the same, since the initial 
field values were essentially exact, but it is reassuring that run C has converged as well toward the same solution. 
Indeed, results for configurations at later times continue to converge toward the exact solution. In Fig. 0] we show 
the evolution of the maximum central density for the three runs, as well as the predicted value from the OV solution. 
The maximum density from run A stays near this value throughout the evolution, with small deviations which result 
from the unavoidable discretization effects present in SPH; in general, SPH will yield very accurate global integrals 
over a mass distribution, since numerical noise smooths out, but demonstrates significant noise in quantities defined 
for individual particles, which vary iteration to iteration as each particle's neighbor list adapts to current conditions. 
The maximum density for run B oscillates around the proper value with a period of T/Mq = 112, showing no signs 
of systematic drift. This is very close to the proper value for the limiting case of infinitesimal radial variations, 
T/Mo = 104, which we find by interpolating from the values given in Table A18 of J65|, after scaling their results to 
our units. 



B. Spherical Dust Cloud Collapse 

To further test the dynamical aspects of the code, we computed the evolution of a uniform density dust cloud, 
i.e., a spherical distribution of matter with zero pressure, started initially from rest. This is a familiar problem from 
cosmology, and the solution is well-known, but the conversion to CF coordinates and the matching conditions at the 
surface of the dust cloud make the resulting expressions considerably more complicated. The complete description 
of the metric as a function of time was derived for the case of both maximal time-slicing |66| and polar time-slicing 
|67| : it is the former case which can be compared to our results. In this approach, the field values and CF time and 
position variable values are computed by solving a set of ordinary differential equations which describe the behavior 
of the fields in terms of the comoving time and space variables. Using these, it is a simple process of interpolation 
to derive the fields as functions of CF time and position. We checked our integration code by comparing our results 
against the plots in [66|, finding perfect agreement. 

To construct our initial configuration, we used essentially the same techniques described above. Particles were laid 
out in an HCP lattice, with masses set proportional to the analytically known value of p*(r), which we derive as 
follows. Note that varies with radius; it is p(r) that is initially uniform. 

The total mass of the cloud was set to unity, and the initial radius in comoving coordinates to r s , the parameter 
used for all figures in |66| . The initial metric in comoving coordinates, for a cloud with unit mass and comoving radius 
f s , is given by 

ds 2 = -dr 2 + a(T) 2 (d X 2 +sm 2 X dn 2 ), (83) 
o(r = 0) = 2/(sin Xs ) 3 , (84) 
f. = 2/(sin Xs ) 2 , (85) 

with the exterior metric described by the Schwarzchild form, Eq. 1721 with f = a(r = 0) sinx. In terms of Xs, the CF 
radius is given by 

r s =l(r s - M + Vr s (r s ~ 2M )) = \ I — 2 l + J — [ r^— j 
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2 1-cos Xs V ' 
Equating interior metric coefficients, we set Adr = ad\ and Ar = asin X , dividing and then integrating to find 

/i \ 1/2 

r oc f^2L] , (87) 

Vi + cos X y 

i/i \ 1/2 / i , \ 3/2 

1 / 1 — cos x \ / 1 + COS Xs 



2 \ 1 + cos xj \ 1 — cos x 

where the proportionality constant is determined from Eg. 1861 The initial conformal factor can now be written down, 
since 

A = asinx = 40 + cosx) (gg) 
r (1 + cosXs) 

and we can simplify the resulting expression by noting that cosx s = (2r s — l)/(2r s + 1) and cos X = (1 — r2 /2 r s)/(l + 
r 2 /2r^). Since the matter starts from rest, we know that Ui = initially, and thus 7„ = 1 everywhere, and we find 
the initial rest-mass density profile is given as a function of radius by 
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Once the initial particle configuration was set, we calculated the dynamical evolution of the system using the same 
techniques described above for the OV case. Since the pressureless material had no outwardly directed force, the 
inevitable fate of the system was collapse to a BH. 

To allow for direct comparison with the figures shown in |66| , we computed the evolution of a dust cloud with unit 
ADM mass and an initial areal radius f s = 10, just as they did. In Fig. \E\ we show the evolution of a set of equally 
spaced Lagrangian tracer particles compared to the exact semi-analytic solution we computed. This corresponds 
to their Figs. 8-9, with two slight differences. Our figures are plotted in CF coordinates, rather than comoving 
coordinates, and our tracers are equally-spaced in radius, not in increments of enclosed mass. For comparison with 
their Fig. 9, wc also show the more familiar areal radius of the cloud, which roughly satisfies the relation r s f s — 1 
initially. We see the agreement between our calculation and the exact solution is very good throughout the evolution, 
up until t/Mo w 38, where a slight discrepancy begins to develop, primarily at the surface of the cloud. This time 
corresponds closely with the formation of an event horizon for the system, shown as a dotted line, starting at the 
center at t/Mo — 38 and moving outward, reaching the surface of the cloud at t/Mo — 43.4, shown as a horizontal 
line. This late time discrepancy has two sources. The first is that SPH, which by definition produces a differentiable 
density field, cannot reproduce the step-function density drop at the surface of the dust cloud. This explains to a 
large degree why the outermost tracers diverge furthest from their exact path. In addition, the large field values and 
gradients found around the event horizon present a challenge for our field solver. We typically see a dramatic increase 
during this period in the number of relaxation iterations required to converge to a sufficiently accurate solution. 

In Fig. El we show the evolution of the conformal factor A and lapse function N for the dust cloud, again in 
comparison to the exact solution, at times t/Mo = 0, 20, 30, and 40. We see again that our code can reproduce 
the proper solution past t/M = 30, but by t/Mo = 40 shows non-trivial deviations from the proper solution. At 
t/Mo = 40, the relative error in the metric fields and position of the Lagrange radii is approximately 4%, growing to 
roughly 15% by the time the event horizon encompasses the entire mass distribution. In all cases, the deviations from 
the correct solution take the same form: our computed field values are closer to unity (and the previous timestep's 
solution) than we would expect from the semi-analytic solution. 

This behavior was confirmed by testing the collapse of dust clouds with initial areal radii of f s = 5 and f s = 50. 
In both cases, we find extremely accurate results until the event horizon forms, at which point our accuracy degrades 
to a noticeable extent. The effect seems to depend primarily on the formation of an event horizon in the system, and 
not on the relative change of the density or various field quantities. We have concluded that the relaxation techniques 
used in the code have difficulty in their current form in handling the steep spatial gradients in the shift function near 
the event horizon (see Fig. 2 of 66]), and we are working on techniques to better handle this situation. Of particular 
importance is altering the relaxation parameters of the iterative scheme used by the field solver in the presence of 
these large field values and gradients near the event horizon, to correct for the systematic drift away from the expected 
values. 
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C. Circular Orbits of Quasi-equilibrium Configurations 

Finally, to test out the overall stability of the code, we evolved several quasi-equilibrium binary configurations using 
our code, but without adding in the radiation reaction drag terms. Since the CF formalism is time-symmetric, we 
expect that any stable quasi-equilibrium model should yield a circular orbit, maintaining a constant binary separation, 
ADM mass, and internal rotation profile, among other parameters. Here, we chose models at an initial binary 
separation of r^j M. c h = 19.91, 20.42, and 22.97 from the M/R = 0.14, 0.14 equal-mass sequence of TG, denoting the 
resulting calculations as runs QC1, QC2, and QC3, respectively. For NS with ADM masses of Mo = 1.4 M Q , these 
separations correspond to 33.64, 34.51, and 38.81 km, respectively. The innermost of these represented the limiting 
case found for this sequence just before formation of a cusp at the inner edge of the two NS. 

As a first check of our code, we compared the orbital frequency determined from our dynamical runs to the known 
value for each model taken from the quasi-equilibrium sequence. We found excellent agreement between the orbital 
periods computed from our runs, T/M ch = 422.6, 438.2, and 506.4, and those determined by TG, T/M c h = 422.8, 
438.0, and 519.7, for runs QC1, QC2, and QC3, respectively. 

There are several conserved quantities which should be respected in the time-symmetric CF formalism, allowing for 
further code tests. In Fig. [7| we show the evolution of the binary separation for the runs. For each run, the orbital 
period is shown with tick marks. We see in each case that the orbit is stable, with variations in the binary separation 
of no more than 4% during the first two orbits. The timescale for the radial variations is similar to the orbital period, 
but not an exact match; this reflects both the effects of GR as well the slight degree of time-asymmetry present in 
the numerical implementation of the CF formalism. We note that the deviations from circularity were largest for 
run QC3, which had the largest binary separation. We believe this results from the larger relative magnitude of the 
spurious initial velocities that result from deviations away from equilibrium in the initial condition; while these terms 
are of essentially constant magnitude in all three runs, the equilibrium velocity field has the smallest magnitude at 
the largest separation. 

Run QC1 was started from the innermost point along this binary NS equilibrium sequence, and the binary performs 
three complete orbits with no sign of plunging behavior. AS such, these results can be taken as the first direct proof 
that the entire equilibrium sequence is stable, and suggest that these configurations should be reasonably accurate 
approximations to the true physical state of merging binaries. Further evidence of this claim is presented below in 
Sec. IV Bl when we describe the results of our calculations with radiation reaction effects included. Note that this 
result is not unexpected given the absence of a turning point (minimum of ADM mass and total angular momentum) 
along this irrotational equilibrium sequence. Indeed, while such a turning point along an equilibrium sequence of 
corotating binaries marks the onset of secular instability, a true dynamical instability is usually associated with a 
turning point along an irrotational sequence 55, 68, 69] . 

While this result may appear at first glance to disagree with those of / citeMarro, who find that the ISCO occurs at a 
greater binary separation than the termination point of an equilibrium sequence, we believe that the difference is purely 
semantic. We define an initial configuration to be dynamically stable if in the absence of dissipative radiation reaction 
effects, the circular orbit remains stable (no merger occurs) when evolved forward in time. In |70j| . a configuration 
is described as within the ISCO if in full GR (i.e., including dissipative radiation reaction effects), a binary starts 
merging (the surfaces of the two NS come into contact) within one orbit when evolved forward in time. Clearly, 
using these definitions, the same initial configuration can be dynamically stable (in the time-symmetric CF sense) at 
a separation within the "ISCO" as defined by [7(j . 

In order to estimate how well our code respects conserved quantities, we show in Figs. |H1 and ED the evolution of 
both the ADM mass and system angular momentum, calculated from the SPH expressions, Eas.1421 and 1441 and the 
spectral basis forms, Eas.l41land l43l respectively. In the former, we see that the SPH expression for the ADM mass is 
remarkably constant over time, with only very minor deviations of relative magnitude less than a tenth of a percent. 
The system angular momentum varies more, but is still conserved to within 1%, with no sign of a systematic drift 
in either direction. We note that much of the variation is correlated with the deviations in the binary separation, 
oscillating on the orbital timescale. We see much more variation on an iteration to iteration basis when we look at the 
same quantities computed using the spectral basis. This is hardly surprising, since these values are computed at the 
end of a relaxation routine, and we expect some degree of variation on a step-by-step basis depending on the exact 
phase space path traced out by the iterative solution. All the same, we see that the angular momentum is conserved 
to roughly the same level in the spectral basis as it was when computed using SPH, and while the variation in the 
ADM mass is larger, we see no sign of a systematic drift. 
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V. DYNAMICAL CALCULATIONS 

While dynamical calculations including the effects of radiation reaction are the only way study the coalescence of 
binary NS systems, they have several additional uses which are often overlooked. Of particular importance is the 
ability to determine the validity of quasi-equilibrium models as initial conditions for dynamical calculations, regardless 
of whether CF or full GR gravity is used. Thus, we computed two dynamical runs including radiation reaction. Run 
RR1 was started from the cusp point of the sequence (the same initial configuration used in run QC1), and was 
used to study the details of the coalescence process. Run RR2 was started from a larger separation (the same initial 
configuration used in Run QC3), and was used to test out the deviations we expect from quasi-equilibrium prior to 
reaching the termination point along the equilibrium sequence. 

A. Stable Regime 

The evolution of the binary separation for run RR2 is shown in Fig. [TJ3 The dotted curve, showing the result 
from the calculation, does not have the behavior one would expect. Notably, after the binary separation decreases 
monotonically until reaching r / Ai c h — 21.2 at t/A4 c h — 700, the system turns around and expands briefly back to 
a separation of r/M. c h = 21.5 at t/JA c h = 950 before shrinking again. This does not reflect any inherent problem 
in our radiation reaction formalism, since oscillations in the binary orbit were found for calculations which ignored 
all radiation reaction effects. In fact, if we "correct" the binary separation by looking at the difference in separation 
at equivalent times between runs RR2 and QC3, which were started from the same initial configuration, with and 
without radiation backreaction terms, we see a pattern of monotonic decrease, shown as a solid line. This seems 
to indicate that deviations from circularity in the orbit of the quasi-equilibrium binary configurations represent a 
systematic effect in the evolution. 

The "corrected" infall curve shows clear signs of an orbital eccentricity with a timescale roughly corresponding to 
the orbital period. This is a natural consequence of starting out from an initial condition with zero infall velocity, and 
has been seen before in virtually every PN and CF calculation (FR, j^liE])- Its origins are clear: the framework used 
by GGTMB to construct quasi-equilibrium initial conditions assumes a helical Killing vector exists, which enforces 
an initial circularity in the orbit, rather than the proper infalling trajectory. If calculations could be started from 
sufficiently large separations, GW emission would cause the orbit to circularize, but the process works slowly, and 
breaks down as the binary makes the transition toward a dynamical merger. It is fair to say that no computational 
scheme can currently be trusted to remain stable over the period required to circularize the orbit. 

B. Coalescence 

We have also computed the full dynamical evolution of a binary system started from the innermost point along the 
quasiequilibrium sequence, denoted run RR1. In Fig. 1111 we show density contours from the system during the binary 
phase, which lasted until t/M. c h — 883. The contours are equally spaced logarithmically, two per decade, ranging from 
density values of M 2 ch p* — 10~ 6,5 — 10 _1 . We recognize a familiar pattern from past PN and relativistic calculations 
(FR, [33, including the development of tidal lag angles as the timescale on which the gravitational field 

evolves becomes comparable to the dynamical timescale. This gives rise to an "off-axis" collision, as matter from the 
inner portion of each star runs along the trailing edge of the other, forming a turbulent vortex sheet (see FR3). Some 
fraction of the mass in these flows eventually crosses through the outer Lagrange point on the opposite side of the 
binary, forming the very low-density spiral arm structures seen at t/A4 c h = 850. In contrast to Newtonian binaries, in 
which angular momentum transfer outward leads to massive spiral arm formation, the very low-density arms formed 
here have velocities much smaller than the escape velocity, and the vast majority of the mass remains gravitationally 
bound to the system. 

At t/ M. c h — 883, shortly before the binary field solver fails to converge, we take our matter and field configurations 
and transform them into the single-star description described in Sec. lIIIBl In Fig.^J we show the particle configura- 
tion at t/A4 c h = 883, which can be described as a bar with two low-density arms trailing off the edges. While it may 
seem inappropriate at first to describe the configuration as an ellipse, we note that the low-density contours, shown 
as dashed lines, do form a much more elliptical pattern than one might at first expect. In interpreting SPH particle 
plots, it is important to remember that low-density particles have large smoothing lengths, implying that the matter 
distribution extends well beyond the apparent sharp edge. The boundaries of the innermost computational domains 
are shown in the figure as heavy solid lines, for both binary and single-star configurations. We see that they align 
rather well, failing to overlap only in low-density regions near the boundary. To confirm the validity of the switch, we 
compare the field solutions for the particles before and after the transition. In Fig. 1131 we show the relative change in 
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the lapse function N (top panel) and conformal factor A (bottom panel), as a function of the x-coordinate. We see 
that in all cases the relative error is < 1%. 

The evolution of the single-body configuration is shown in Fig. ^] We see a strong pattern of differential rotation 
in the merger remnant, which slowly relaxes from a very elliptical shape toward a more spheroidal one. This is to 
be expected, as our choice of EOS with T = 2.0 should not be able to support a long term ellipsoidal deformation 
[Hsj . By t/Mch — 1220, the remnant has relaxed to a nearly circular profile, but the differential rotation, as shown 
in Fig. 1151 persists. The latter will dissipate sl owly on either the viscous or magnetic braking timescales, beyond the 
scope of what we can reasonably calculate |7l|. l72t l73| . Differential rotation is expected to stabilize the star against 
gravitational collapse in the short-term |7J, [75| . Quantitatively accurate determinations of the rotational velocity 
profile in terms of the parameters of the initial system are likely to be crucial for making prediction as to which 
systems will or will not collapse promptly to BHs, especially since these systems will likely have very large masses. As 
we see in Fig. 1161 the vast majority of the rest mass of the system ends up in the merger remnant itself, with a fraction 
of a percent of the total mass forming a low-density, bound halo around the remnant. This seems to be the consensus 
from PN and relativistic calculations of irrotational binaries (FR3,0|) and even PN and relativistic calculations of 
synchronized binaries (FR3, 48]), which traditionally yielded significantly higher mass ejection fractions in Newtonian 
calculations. We note that PN calculations of irrotational binaries yielded an ejected mass fraction of < 1% (FR3), 
rather than the 6% quoted by |43j . 

To compare with the relativistic results of j we show the evolution of the maximum density as a function 
of time in Fig. 1171 We see at the earliest times a low-amplitude pulsation, resulting from small deviations from 
equilibrium in the initial SPH particle configuration. This pulsation damps away almost completely by the time the 
stars plunge inward. As the system begins to accelerate rapidly inward prior to the merger itself, the maximum 
density decreases as the stars are tidally stretched. This effect, seen in a number of calculations, further indicates 
that these systems will not undergo a pre-merger gravitational collapse. This process was originally suggested in |45| . 
but is likely to have been a numerical artifact since they used a version of the CF formalism containing an error in 
one of the evolution equations. Indeed, the quasi-equilibrium sequences in TG show a slow decrease in the central 
baryonic and energy densities, /?» and p, as the binary separation decreases. After contact, we see a str ong rise in the 
maximum density, followed by a rapid, high-amplitude oscillation. This result is similar to that seen by |48j . although 
they found a relatively higher maximum density value during both the peak and the trough of the oscillation. This is 
almost certainly a result of using different initial spin configurations. Irrotational systems, such as the one used in run 
RR1, concentrate relatively more angular momentum at small radii compared to initially synchronized systems, like 
run A of |4S[ . Thus, when a remnant is formed from an initially irrotational binary, the central density will typically 
be lower, since there is a greater centrifugal barrier and less pressure support is needed to stabilize the configuration. 
Our results differ rather significantly from those of |43|. who found a smaller-amplitude, more sinusoidal oscillation 
after merger. It is possible that this discrepancy can be attributed to the use of the CF approximation rather than 
full GR. A much more likely explanation, however, is that the difference results from the numerical methods used 
to describe shock heating in the matter. Lagrangian SPH codes were used here and in [lif. whereas Shibata and 
collaborators 0, 0> ^1 use an Eulerian grid-based code, which may have a better ability to resolve shock fronts. 
It seems clear by examining the results from run M1414 in [43 that as the NS cores converge, the increase in the 
central density is suppressed by the conversion of kinetic energy into heat. The merger remnant shows a spike in the 
internal energy in the very center, 2 — 3 times the adiabatic value, and a corresponding decrease in the density, giving 
the remnant its toroidal profile. Only calculations using their new shock-capturing scheme |76j produce this behavior; 
previous calculations yielded remnants with centrally condensed density profiles 0, • ^ will be of great future 
interest to determine whether or not SPH studies of rapidly rotating collapsing matter configurations can produce 
these toroidal configurations with "hot" cores, and if so, which aspects of shock physics are crucial for understanding 
this process. 

In Fig.^0 we show the GW signal for run RR1, calculated from Eas. 1551 and 1561 The waveforms show a familiar 
chirp signal up until t/M. c h ~ 850, followed by a modulated, high-frequency ringdown component. The strength of 
our signal at peak amplitude matches extremely well with the results of [43 . l48j |. as one would expect from simple 
dimensional analysis. Our modulated remnant signal, though, is much more similar to the results of |43j| than 
|48|. who find a damped ringdown signal of lower amplitude for this model, with no obvious modulation. Previous 
PN calculations (FR2, FR3) identified the source of the GW amplitude modulation as a combination of differential 
rotation and ellipticity in the remnant. When the inner and outer regions of the remnant have ellipsoidal deformations 
which are roughly aligned, as we see at t/A4 c h ~ 900 in Fig. ^2 the GW amplitude will be at a maximum. When 
differential rotation drives the inner regions into misalignment, as we see at t/M. c h sa 1020 in the same figure, the 
GW amplitude reaches a minimum. Eventually, these effects dissipate away as the remnant relaxes toward a more 
spheroidal configuration. We believe the lack of modulation in the waveforms shown in Fig. 10 of JZ3 merely reflects 
the use of a different initial spin configuration. As stated above, synchronized initial configurations contain relatively 
less angular momentum at smaller radii, and the lack of a centrifugal barrier allows the remnant ellipticity to damp 
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away much more quickly as the densest regions in the NS cores fall into the remnant's center. Our results are in broad 
agreement with those of |43j , who find a modulated waveform immediately after the merger on a similar timescale as 
our results, before a longer-lived, smaller-amplitude damped modulation eventually appears. 

Calculating the GW energy spectrum for the merger waveform is more complicated than simply applying Eq. 1581 
since taking the Fourier transform of a signal with non-zero initial and final values introduces aliasing of the boundary 
conditions into the resulting waveform. To correctly derive the proper spectrum, one must also "attach" analytic 
solutions to the beginning and end of the calculated signal, representing the portions of the inspiral and ringdown 
phases, respectively, which fall outside the bounds of the numerical evolution. In the past, the authors (FR3), and 
others [Ig, |63, |Z3 > have modeled the initial inspiral phase via the Newtonian point-mass approximation, but that is 
clearly not physically realistic. Indeed, it has been shown in FGRT that relativistic effects should have a significant 
effect on the GW spectrum even before the dynamical merger. Summarized, we know that for an equal-mass binary 
system, the total mass-energy is given by E = M t — M 2 /8r 2 — 2 12 M c h — M 2 c hl '(2 a6 r 2 ), and the (Keplerian) orbital 
frequency by fxe P = (\/Mt/r 3 )/2Tt. In terms of the GW emission frequency, few = 2/^ep = ("2 a ' 6 /n)y/-M c h/r 3 , we 
find 
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where the latter equation demonstrates the familiar power-law dependence of the GW energy spectrum. 

For the quasi-equilibrium sequence from which we take our initial condition, we found that the system ADM mass 
can be given in terms of the GW frequency by a (phcnomcnological) fit of the form 

E(f G w)/M ch = E N /M ch - 0.4905(.M C ,JW) + 23l(M ch f G w) 2 ■ (93) 

Differentiating this equation yields what we term the "quasi-equilibrium energy spectrum" , but we require additional 
assumptions to be made before we can construct the time-history of the inspiral waveform. First, we determine a fit 
for the GW frequency as a function of conformal separation, finding that we can approximate the proper function to 
within 0.1% with the form 

Mchfaw = —(r/Moh)- 1 - 5 - 1.592(r /M ch y 2 - 5 + 6.325(r/M ch )- 3 5 . (94) 
Next, we assume that the GW signal amplitude is given by a slightly modified version of the quadrupole form, 

QiW = ~Q l yl(t) = (1.0 - K )2°- 2 TT 2 M ch r 2 f GW co S (9 G w(t)), (95) 
Q!|(i) = (1.0-K)2 - 2 ir 2 M ch r 2 f GW sin(9 G w(t)), (96) 

where the prefactor k is used to match the amplitude of the inspiral signal onto that at the beginning of our calculated 
signal. We find that n — 0.015 throughout the early phases of our calculated waveform, indicating that it can be 
well-approximated by the expected quadrupole form. Furthermore, we note that the resulting energy spectrum is 
essentially independent of k; as k increases, the GW amplitude decreases, decreasing both the energy loss rate and 
the frequency sweep rate in the same proportion, leaving dE/df unchanged. Finally, we assume that the energy loss 
rate in GWs is given by the standard quadrupole expression, 

dE/dt = 0.2 (4n 2 f GW Q [ §Q [ §.) (97) 

To construct the inspiral waveform, we start from a point along our calculated waveform and evolve backward in time 
from that point. At every timestep, we calculate the instantaneous energy loss rate from Eq. E| After adjusting the 
total energy, we calculate the new GW frequency by implicitly solving Eq.|^l and adjust the phase of the GW signal 
appropriately. We find the new binary separation by solving Eq. 1941 implicitly, and finally evaluate the waveform via 
Eqs.l5]andOS] 

The question of where to match the quasi-equilibrium waveform to the calculated one deserves some attention. 
Matching the two at t — is very much a mistake, because it represents a transition from an infalling configuration 
to a circular one; beforehand, the frequency sweep rate dfew/dt is positive and increasing, while afterward it is reset 
instantaneously to zero. This mismatch in the infall velocity, and thus the frequency sweep rate, results in energy 
"piling up" at the transition frequency, as can be seen in FR3 and to a smaller degree in Fig. 12 of |48f. We find that 
matching the inspiral waveform to the calculated signal at t = reproduces this error, but that by t/M. c h — 250, the 
match in the inspiral velocity is sufficient to leave no measurable trace in the resulting spectrum. 
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We note that |48( also match their inspiral waveform to their calculated one at some time into the calculation, 
but we believe that they place too much trust in the behavior of the energy spectrum near this transition frequency. 
In their paper, they alter the frequency of a Newtonian inspiral waveform to match their relativistic calculation by 
adjusting by hand the coalescence time, r = (dr/dt)/r. We believe this approach to be a mistake. The effect of 
this change is to use a Newtonian waveform with different physical parameters from those used in the calculation; in 
particular, it is equivalent to calculating a Newtonian waveform using the wrong total mass, and thus the wrong limit 
for the spectrum at low frequencies. It is only through the use of an inspiral waveform whose frequency approaches the 
proper Keplerian limit at low frequencies and the relativistic form at high frequencies that accurate energy spectra 
can be constructed. Indeed, since our initial condition was taken from the same quasi-equilibrium sequence used 
to generate our frequency data, the initial GW frequency derived from our calculation matched that of the inspiral 
waveform within 0.1% without requiring further manipulation. 

In Fig. ^3 we show as a solid line what we believe to be the first complete and consistent relativistic waveform for 
a binary NS merger at frequencies few < 1.5 kHz. The frequencies listed on the upper axis assume our "standard 
model" parameters, i.e., each NS has an ADM mass Mq — IAMq. The two dotted lines show the components 
which make up the energy spectrum. At low frequencies, Mchfcw < 0.004, we see the primary contribution is 
from the quasi-equilibrium inspiral waveform, whereas at higher frequencies it is from our calculated waveform. The 
short-dashed line shown the Newtonian point-mass relation, given by Eq. 1921 and the long-dashed curve the quasi- 
equilibrium result, found by numerically differentiating Eq. 1931 We see excellent agreement between our calculated 
waveform and the quasi-equilibrium fit, up until frequencies M c hfcw ~ 0.007 — 0.009. This peak represents the 
"piling up" of energy at the frequency corresponding to the phase of maximum GW luminosity, as the stars make 
contact and the infall rate drops dramatically. The second peak, at Mchfcw ~ 0.010 — 0.011, represents emission 
from the ringdown of the merger remnant. It is likely that we underestimate the true height of this second peak, since 
we assume that the GW signal after our calculation damps away exponentially. Still, it is extremely unlikely that 
including the ringdown phase will increase the strength of this peak by more than a factor of a few, since our chosen 
EOS will not support a long-lived ellipsoidal deformation, and it is likely that the ringdown oscillations will smear 
the GW emission over a small range of frequencies rather than coherently emitting at a single frequency. 

In general, the energy spectrum we calculated here confirms the general conclusions we put forward in FR3 and 
FGRT, albeit in a much more consistent way. The GW energy spectrum does show a significant drop away from the 
Newtonian point-mass form at frequencies significantly below 1 kHz, in almost the exact same form as we predicted 
from quasi-equilibrium data alone in FGRT. Nowhere does the spectrum rise above the Newtonian value, including 
the peaks associated with maximum GW luminosity and the ringdown oscillations. These results suggest that the 
weak signal amplitude of the peaks above 1 kHz, which lie outside the Advanced LIGO broadband frequency range, 
may inhibit detections by high-frequency narrow-band interferometers as well. However, combining lower-frequency 
narrow-band detectors with broadband LIGO measurements, as suggested in |34|, appears extremely feasible, and 
may allow GW measurements to constrain the NS compactness and EOS. 

One possible cause for concern with our code is non-conservative behavior caused by numerical errors that develop 
after t/A4 c h = 500. In Fig. [201 we show the change in the system's angular momentum over time. The dotted curve 
is the uncorrected result derived from our calculation, which shows two periods of angular momentum generation, the 
first from t/A4 c h = 600 — 700, and the second at t/A4 c h = 875 — 900. The former is associated with the numerical 
inaccuracies discussed in Sec. IV Al for run RR2. All of our runs show some spurious angular momentum generation 
and slowing of the binary infall during this period. The latter spike occurs immediately before and after the transition 
from a binary to a single-star description. Correcting for both of these spurious terms yields the solid line on the plot, 
which still underestimates by a non-trivial amount the angular momentum loss we would expect from the quadrupole 
formula, 

{j k )GW = QAe l3k (Q$Q [ ^). (98) 

Using the qua drup ole formula on our results yields a total angular momentum loss fraction which very nearly equals 
that found by |43,[48|]> with approximately 7% of the system's angular momentum converted into GW emission. The 
discrepancy between this amount and what we derive from the SPH particle configuration, Eq. 1441 can be easily 
understood. First, we do not see strong angular momentum loss as the NS first make contact, since this is where 
we push our field solver to its limit. Second, our method for estimating the instantaneous angular velocity, Eq. 1541 
underestimates the proper value of w, yielding a radiation reaction force, Eq,021 smaller than the correct quadrupole 
value (which we can determine after the fact). This error could be decreased in magnitude by calculating the angular 
velocity from the change of position of the NS centers-of-mass in time, but such a prescription is difficult to define 
consistently in the single-body regime. Even defining the orbital frequency in terms of the rate of change of the 
quadrupole tensor, as done by |4S| , underestimates the correct angular momentum loss rate by up to 40% during the 
GW emission peak. 
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We see a similar pattern at work in the evolution of the system's ADM mass, shown in Fig. comparing the 
energy loss to GWs from the quadrupole approximation formula (dashed line) , Eq. %J7\ to the value we find from SPH 
summation via Eq. 021 (solid line). The quadrupole value agrees well with other calculations, which typically find 
AM adm /M adm = 0.004. Looking at our particle summation value, we see a slow decrease from t/M. c h = — 500, 
of magnitude 0.2 — 0.3% of total value, in good agreement with other calculations and our own quadrupole estimate. 
From that point on, we see a small spurious increase from the point the numerical errors begin to become significant, 
followed by a sharp decrease of ~ 2% immediately prior to our transition from binary to single-star descriptions. 
Once we have made the transition, the ADM mass oscillates slightly, with an overall peak-to-trough amplitude of 1%. 
Thus, we conclude that while the instantaneous value we measure from the particles directly is liable to be off by up 
to 2%, we can reconstruct the proper energy loss rate after the calculation is over. 

VI. CONCLUSIONS 

We have developed and tested a new relativistic 3-d Lagrangian hydrodynamics code, which should prove useful 
for studying a wide variety of physical systems. Here, as an initial investigation, we have performed the first full 
evolutions of the coalescence and merger of irrotational binary NS in the CF approximation to GR. Moreover, these 
calculations represent the first numerical evolution of coalescing binary systems performed with either a spectral 
methods field solver or the use of spherical coordinates adapted to a binary environment. 

The code has been validated using several tests. We can accurately reproduce static spherical stellar configurations 
as well as the known solution for a collapsing pressureless dust cloud. In both cases, the CF approximation yields 
the exact solution in GR, as it will for any spherical configuration. We find that we can reproduce these known 
semi-analytic solutions to high accuracy, up until the formation of an event horizon for the collapsing dust clouds. 

Our dynamical evolution calculations for quasi-equilibrium models at a number of different binary separations 
indicates that we can successfully integrate forward for several orbits, with typical errors in conserved quantities of 
~ 1%. In doing so, we have demonstrated directly for the first time that the M/R = 0.14, equal-mass sequence of 
TG is stable all the way to its innermost configuration, at which point a cusp develops on the inner edge of each NS. 

Our dynamical calculation of a complete binary NS merger, including radiation reaction effects, demonstrated that 
our spherical coordinates, spectral method approach is robust enough to follow the system from the point just before 
the formation of a cusp through merger and the formation of a stable remnant. Some errors were introduced during 
this period into the globally conserved quantities such as the ADM mass and system angular momentum, but we 
find that the field values were computed consistently throughout, and that the global dynamics was treated in a 
quantitatively accurate way. 

We find that the merger remnant formed in our calculation is differentially rotating, with a transient quadrupole 
deformation. This combination of effects produces a GW amplitude with a modulated form, similar to what has 
been seen before in PN calculations (FR3) and more recent full GR calculations of the same model |43|. We find 
that the remnant is initially stable against gravitational collapse, as did |43|. with the supermassive NS (which has a 
baryonic mass essentially twice that of either NS in isolation) supported by strong differential rotation. We find that 
a density maximum develops rather rapidly in the center of the merger remnant, as has been seen in all other PN and 
CF calculations, but not that of 0], whose full GR merger calculation yielded a toroidal remnant. We believe the 
difference results from their use of a capturing scheme, whereas our runs were performed using an adiabatic treatment. 

By combining our calculated GW signal with a relativistic quasi-equilibrium inspiral precursor, we have generated 
the first GW energy spectrum from a binary NS merger which is complete at all sub-kHz frequencies and consistent 
throughout. We find that the energy spectrum deviates from the Newtonian power-law relation by more than 50% at 
frequencies few < 1 kHz (the "break frequency"), in very good agreement with the predictions of FGRT. There are 
distinct peaks in the power spectrum corresponding to the phases of maximum GW luminosity and merger remnant 
ringdown, but at levels significantly below the point-mass power-law value. 

In our future work on binary NS systems, we hope to address a number of topics, many of which deserve much 
more careful study. Based on the excellent agreement between our calculated GW energy spectrum and that based 
purely on equilibrium sequence data, we hope to do a broad phase space survey to determine the dependence of the 
"break frequency" on both the NS EOS and the system's mass ratio. Beyond this parameter study of NS-NS mergers, 
we also plan to investigate in detail the formation process for the merger remnant, to determine the conditions which 
may lead to the formation of a quasi-toroidal merger remnant. This will necessarily involve the use of a relativistic 
artificial viscosity scheme to treat shocks. The density profile of the merger remnant is likely to influence the final fate 
of the system, and may prove crucial for determining the coincidence properties of GW emissions and short-period 
GRBs, should they result from compact object binary mergers, as has been widely suggested. 
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TABLE I: A comparison of our notation for various relativistic quantities to previous works using the CF formalism: prl liH 
I47I. Efij. For those cases where no unique terminology was defined, we give the simplest equivalent algebraic form. We also list 
the equation in this paper where the quantity is defined 



Quantity 


Here 


Gourgoulhon 


Oechslin 


Shibata 


Wilson 


See Equation 


Lapse 


N 


N 


a 


Q 


a 


m 


Shift 


Ni 


Ni 


-ft 


-ft 


-ft 


m 


Conf. Fact. 


A 


A 


V> 2 


v> 2 


2 


m 


Rest dens. 


p, 


r n A 3 P 


p* 


p* 


Dff 


□ 


Lorentz Fact. 


Jn 




au° 




W 


m 


Velocity 


v* 


NU* + N l 


v i 


V* 


V 1 


m 


Spec. Momentum 


Ui 


Wi 


Ui 


Ui 


Si/(D4, 6 ) 


ma 


Enthalpy 


h 


h 


w 


1 + T6 


h 


mi 



TABLE II: A summary of the runs presented in this paper. 



Run 


Description 


See Sec. # 


OVl 
OV2 
OV3 


Equilibrium OV model, GMo/Rc 2 = 0.126 
Non-equilibrium OV model, Ro = l.lR eq 
Same as OV2, w/relaxation drag 


HvAl 

|iVA| 
LvAl 


DC5 
DC10 
DC50 


Collapsing dust cloud, Mo = 1, Ro = 5 
Collapsing dust cloud, Mo = 1, Ro = 10 
Collapsing dust cloud, Mo = 1, Ro = 50 


|1YB| 
|IVBl 

hvfh 


QCl 
QC2 
QC3 


Quasi-circular binary orbit w/o Rad. Reac, ro/M c h = 19.91 
Quasi-circular binary orbit w/o Rad. Reac, ro/M c h = 20.42 
Quasi-circular binary orbit w/o Rad. Reac, ro/M c h = 22.98 


pvcl 

nvn 


RRl 
RR2 


Full binary evolution w/rad. reac, ro/Mch = 19.91 
Full binary evolution w/rad. reac, ro/Mch = 22.98 


IVBl 

EH 




FIG. 1: Radial domains used to evaluate the field equations of the CF method. The boundaries of the inner two domains are 
shown as lattices, with all collocation points within these domains shown as well. The outermost domain, which extends to 
spatial infinity, is not shown. 
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FIG. 2: A pictorial demonstration of the coordinate transformations described in Secs. lIIl'Al and lni Bl the particle configuration 
is for demonstration purposes only, and was not taken from our calculations. In the upper left panel, we show a NS-NS binary 
in the x — y (inertial) frame, as well as the x' — y frame, defined so that the centers-of-mass of the stars (crosses) lie equidistant 
from the origin on the x'-axis (see Eqs. I25H29H . In the upper right, we see the same system in the x' — y frame. The angles 
$i and $2 are determined from the respective moment of inertia tensors using Eq. 1301 The best fit ellipsoidal configurations, 
determined from Eqs. I31H38I are shown around each star, aligning with the X' q — Y' axes, determined from Eqs. IbOl and 1611 
In the bottom left, we show isocontours for the radial functions /i and /2, defined by Eo. 1591 as well as the boundary of the 
overlap region (heavy solid line). Finally, in the bottom right, we show the rescaling of the surface function for very close 
configurations, showing only star 2 for clarity. Here, the binary separation is r/M ch — 6.0, implying the maximum extent of 
the surface of star 2 is to x' /Ai c h = r/4A4 c h = 1.5. We linearly rescale the surface function, as well as the corresponding 
values of f^, for all points with x > 0. 
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3 4 6 

FIG. 3: Conformal factor A and lapse function N for an Oppenheimer-Volkov solution with a F — 2 poly tropic EOS, and 
conformal radius r s /Mo = 6.874 (solid lines), compared to our computed values at t/Mo = 100 for two models: run OV1 
started from equilibrium (crosses) and run OV3 started from a configuration 10% larger with a velocity damping term used 
to drive the system toward equilibrium (triangles). The agreement is within 1% for all particles. Units are defined such that 
G — c — 1. All masses, lengths, and times in the OV and dust cloud calculations are made dimensionless by scaling results 
against the system's initial gravitational (ADM) mass. Note that the conformal radius is not equivalent to the areal radius 
typically used in solving the OV equation. 
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FIG. 4: Evolution of the maximum density for three runs based on the OV solution described in Fig. Run OV1 (solid 
line) was started from equilibrium, and shows only small variations which result from typical uncertainties in SPH summation. 
Run OV2 (dotted line) used the same EOS, but was started from a radius 10% smaller than the equilibrium value, showing 
sinusoidal oscillations with a period of T/Mo = 112 and no systematic drift. Run OV3 (dashed line) was started from the same 
configuration as run B, but with an overdamped drag term to force the system toward equilibrium, converging rapidly toward 
the proper value of M$p ma x = 0.006. 
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FIG. 5: A comparison between the actual paths traced out by a set of equally-spaced Lagrangian tracers in run DC10, our 
calculation of a collapsing dust cloud with unit ADM mass and initial areal radius r s = 10, (dashed lines) and the exact 
semi-analytical solution of |o6|. shown as solid lines. All radii are shown here in conformal coordinates. We see excellent 
agreement up until the point where an event horizon forms at the center of the cloud, at t/Mo = 38. The event horizon moves 
outward (heavy dotted line), eventually enclosing the entire cloud at t/Mo = 43.4, shown as a horizontal line in the figure. At 
this point, when the matter can be properly defined as a BH, our field solver stops converging. For comparison with Fig. 9 of 
[gfj, we also show, as a long-dashed line, the exact solution for the cloud's surface in comoving (areal) radii, f s (t). 
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FIG. 6: The evolution of the conformal factor A and Lapse function N for the dust cloud in run DC10, compared to the exact 
solution. We see, at t/Mo = 0, 20, 30, and 40, the SPH particle values for the lapse (the four curves with values less than 1.0) 
and conformal factor (values greater than 1.0), shown as points, and the exact solutions, shown as dashed lines. The agreement 
is good up until t/Mo = 30, but at t/Mo = 40 we see some quantitative disagreement, since the field solver breaks down as we 
near the point where the cloud collapses completely into a BH. 
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FIG. 7: Binary separation as a function of time for the binaries in runs QC1, QC2, and QC3, using initial data generated by 
I??!, evolved forward in time with radiation reaction terms neglected. The calculations use binaries with initial separations of 
ro/Mch = 19.91, 20.42, and 22.97, respectively. Full orbital periods, with T/M c h = 422.8, 438.0, and 519.7, respectively, are 
shown with tick marks. We see that the resulting orbits are nearly circular, with changes in separation of no more than 4% 
over the first two orbits. This is the strongest available evidence that the innermost point along this equilibrium sequence, is 
actually stable against merger. 




FIG. 8: ADM mass and total angular momentum, calculated from Eas. 1421 and 1441 runs QC1, QC2, and QC3. We see that the 
ADM mass is conserved almost exactly, and the angular momentum to within 1%, with the variation occurring primarily on 
the orbital timescale. 
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FIG. 9: ADM mass and total angular momentum, calculated in the spectral basis from Eas. 1411 and 1431 for the runs shown in 
Fig. |7| We see roughly the same amount of variation in the angular momentum as was found for SPH summation in Fig. |H| 
The ADM mass shows considerably more variation than the SPH version, but remains well within 1% of the original value with 
no systematic drift. 
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FIG. 10: Binary separation over time for run RR2, started from an initial separation of ro/M c h = 22.97, the same initial 
configuration that was used for run QC3. The dotted line shows the original, "uncorrected" result, including a separation 
increase from t/M c h = 700 — 950. This is primarily due to oscillations associated with numerical noise and deviations from 
equilibrium in the initial configuration. Correcting for deviations from circularity in run QC3, which ignored radiation reaction 
(dashed curve), yields the "corrected" result, shown as a solid line. We see monotonic decrease in the separation over time, 
with an ellipticity induced by our initially circular orbit. 
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FIG. 11: The evolution of the matter run RR1, started from just outside the separation where a cusp develops. We followed 
the evolution through the merger and formation of a remnant. Density contours are logarithmically spaced, two per decade, 
ranging from M 2 ch p* = 1CF 6 ' 5 — 1CT 1 . We see the development of significant tidal lag angles at t/M c h = 500, followed by 
an "off-center" collision. This process leads to the formation of a vortex sheet and a small amount of matter ejection by 
t/Mch = 850 from matter running along the surface of the other NS. 
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FIG. 12: Detailed view ol the binary to single-star transition performed during run RR1 at t/Mch = 883. We show the 
positions of all SPH particles, as well as density contours, shown as dashed lines, logarithmically spaced two per decade. The 
surface of the inner computational domain in both the binary and single-star representations are shown as heavy solid lines, 
with good agreement between the two. 



37 



o 

S-H 

H 

CD 
> 



0) 



0.006 
0.004 
0.002 


0.008 
0.006 
0.004 
0.002 



AN/N 



H 1 h 



AA/A 



J I L 



H 1 1 h 



H 1 1 h 



H 1 h 



J I I L 



J I I L 



J I L 



10 







10 



FIG. 13: The relative change in the lapse function N (top panel) and the conformal factor A immediately preceding and 
following the conversion from binary to single-star representations during run RR1 for every particle, shown as a function of 
the particle's position in the x-direction. The maximum error is approximately 0.8%, with a mean difference of ~ 0.2%. 
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FIG. 14: The evolution of the matter in run RR1, after the transition to a single-star representation, following the same 
conventions Fie, II II We see that a dense remnant forms in the center of the system, surrounded by a thin halo. Some ellipticity 
is seen shortly after the merger, but the system quickly relaxes toward a spheroidal configuration, with maximum density in 
the center of the system, unlike the toroidal configuration found by |43| . 
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FIG. 15: Angular velocity of the merger remnant of run RR1 at t/ Mch = 1220, shown as a function of cylindrical radius, 
r cy i = x' 2 + y 2 . We see strong differential rotation, with the highest angular velocity in the center, decreasing monotonically 
with radius. 




FIG. 16: Enclosed mass as a fraction of the total mass of the merger remnant in run RR1, at t/Ai c h = 1220, expressed as a 
function of (spherical) radius. All but a few percent of the total mass of the system forms the body of the merger remnant, 
with no more than a small fraction of a percent ejected from the system. 
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FIG. 17: Maximum density as a function of time for run RR1. We see the SPH configuration oscillates slightly around 
equilibrium, decreasing slowly as the binary plunges toward merger. During the merger, we see a sharp decrease, followed by 
a large spike upward and evidence for sharp, non-sinusoidal oscillations. 
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FIG. 18: The GW signal in the h+ and hx polarizations for run RR1, as seen by an observer situated along the vertical axis, 
following Eqs. 1551 and 1561 We see a chirp signal followed by a modulated ringdown spike. The modulation is caused by the 
alignment between quadrupole deformations in the inner regions of the remnant core and and those at larger radius. When 
there is strong misalignment, there is destructive interference and the signal amplitude drops, as we see at t/M ch = 920 and 
t/Mch = 1220 in Fig. 1141 . At t/M c h = 1080 the density contours are more aligned, and the amplitude reaches a temporary 
maximum. 
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FIG. 19: GW energy spectrum, M~£dE/df, as a function of the GW frequency, Mchfaw, for run RR1. The dotted lines 
show, respectively at high and low frequencies, the components contributed by our calculated signal and the quasi-equilibrium 
inspiral component. Also shown are the Newtonian point-mass energy spectrum (dE/dfaw^Newt (short-dashed line), Eg. 1921 
and the quasiequilibrium fit (dE/dfcw)QE derived from Eg. 1931 We see confirmation that the "break frequency" calculated 
from a fit of E(f) for the equilibrium sequence (i.e., the frequency at which the energy spectrum decreases to a given fraction 
of the Newtonian level) is reproduced by a full numerical evolution. On the upper axis, we show the corresponding frequencies 
in Hz assuming the NS each have a mass Mo = 1.4M© The two peaks correspond to the phases of maximum GW luminosity 
and ringdown oscillations, respectively. 
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FIG. 20: Angular momentum loss for the binary in run RR1 measured using the SPH integral, Ea. 1441 (dotted line), and the 
result after correcting for spurious angular momentum creation during the binary to single-object transition (solid line). We 
see that the result yields an angular momentum loss approximately half that we would have predicted from the quadrupole 
formula, Eg. 1981 (dashed line), primarily because our estimate for the system's angular velocity used in the backreaction force is 
systematically lower than the GW signal would indicate. The quadrupole result we derive shows that about 7% of the system 
angular momentum is emitted in GWs, in line with previous relativistic estimates. 
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FIG. 21: The evolution of the total ADM mass for the system in run RR1, calculated using the SPH integral, Ea. 1421 f solid 
line). Up until t/A4 c h = 600, we see a slow decrease as energy is emitted in GWs, at approximately the rate predicted by the 
quadrupole formula, Eg. 1971 (dashed line). Beyond t/M c h = 600 numerical inaccuracies lead to oscillations of total amplitude 
~ 2%. From the quadrupole estimate, we find that 0.4% of the system's energy is radiated away as GWs, in line with previous 
estimated from relativistic calculations. 



